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ABSTRACT 

The saturation level of the magneto-rotational instability (MRI) in a strongly radiation dominated 
accretion disk is studied using a new Godunov radiation MHD code in the unstratified shearing box 
approximation. Since vertical gravity is neglected in this work, our focus is on how the MRI saturates 
in the optically thick mid-plane of the disk. We confirm that turbulence generated by the MRI is 
very compressible in the radiation dominated regime, as found by previous calculations using the 
flux-limited diffusion approximation. We also find little difference in the saturation properties in 
calculations that use a larger horizontal domain (up to four times the vertical scale height in the 
radial direction). However, in strongly radiation pressure dominated disks (one in which the radiation 
energy density reaches ~ 1% of the rest mass energy density of the gas), we find Maxwell stress from 
the MRI turbulence is larger than the value produced when radiation pressure is replaced with the 
same amount of gas pressure. At the same time, the ratio between Maxwell stress and Reynolds stress 
is increased by almost a factor of 8 compared with the gas pressure dominated case. We suggest 
that this effect is caused by radiation drag, which acts like bulk viscosity and changes the effective 
magnetic Prandtl number of the fluid. Radiation viscosity significantly exceeds both the microscopic 
plasma viscosity and resistivity, ensuring that radiation dominated systems occupy the high magnetic 
Prandtl number regime. Nevertheless, we find radiative shear viscosity is negligible compared to the 
Maxwell and Reynolds stress in the flow. This may have important implications for the structure of 
radiation dominated accretion disks. 

Subject headings: (magnetohydrodynamics:) MHD — methods: numerical — radiative transfer 



1. INTRODUCTION 

^ ' The inner regions of accretion disks around compact 
CO object are expected to be radiation pres sure dominated 
(N ; (IShakura fc SunvaevlflQTl IPringldflQSl . Near the Ed- 
. dington limit, the radiation field not only dominates 
the energy budget of the accretion disk, but also car- 
' ries a significant fraction of the total momentum. The 
I large accretion rate required to explain the observed 
CO ' X-ray luminosities of some compact objects cannot be 
T—i \ explained by ordinary particle viscosity. Because the 
ILJ , physical p rocess to transfer a i igular momentum was 
• unknown, IShakura fc SunvaevI (|1973f) proposed the a 
^ disk model based on the assumption that the stress 
^ , scales with the vertically integrated total pressure with 
■ a constant coefficient a. Magneto-rotational instabil- 
ity (MRI) is now generally believed to be the physi- 
cal mechanism to produ ce the required stress in ion - 
ized accretion disks I'e.g.. IBalbus fc Hawlevlll991l I1998D . 
For the case with pure gas pressure, numerical simula- 
tions of MRI (e g.,lH awlev et al. 1995; Stone ct al. 1996; 
iMillir fc Stonel[200l show that it saturates as a turbu- 
lent state and the Maxwell stress from the turbulence 
can reach ~ 1% — 10% of the gas pressure. It is there- 
fore important to understand how the MRI saturates in 
a radiation dominated flow. In particular, when the radi- 
ation energy density becomes a non-negligible fraction of 
the rest mass energy density of the gas, the properties of 
the MRI turbulence in the mid-plane of accretion disks 
may be affected by radiation, and radiative viscosity may 
contribute t o the transport of an gular momentum in this 
regime (e.g., ILoeb fc LaoJ[l99l . 



A first step in studying the MRI in radiat i on pr essure 
dominated disks was made bv lTurner etBl\ (|2003| ). who 
studied MRI in the mid-plane of radiation dominated 
accretion disk using a flux-limited diffusion (FLD) mod- 
ule in ZEUS (Turner & Stone 2001) and neglecting ver- 
tical gravity. The vertical structure of a stratified radia- 
tio n dominated a c cretion disk wa.s sub seque ntly studied 
bv iTurned (pOOl . iHirose et al.l (poM ) and iBlaes et al.l 
One limitation of FLD is that it is based on sin- 
gle moment closure, and it drops radiation inertia in the 
momentum equation. When the accretion rate is near 
the Eddington limit, the photon momentum is a signif- 
icant fraction of the fluid momentum, the FLD closure 
may not be appropriate to model the flow. 

Recently, we have developed an impr oved algorithm for 
radi ation MHD (seeljiang et al.ll20ll thereafter JSD12 
and iDavis et al.l I20i2t ). based on a two- moment clo- 
sure that solves the time dependent radiation momen- 
tum equations, and therefore retain radiation inertia. In 
this paper we use this algorithm to study the satura- 
tion of the MRI in the radiation dominated accretion 
disks. As a first step, we focus on the mi d-plane of an 
accretion disk neglecting vertical gravity. [Turner et al.l 
(|2C)03,) found that for the parameters corresponding to 
the radiation dominated regime of the standard a disk 
model, the turbulence generated by MRI becomes highly 
compress ible, consistent with exp ectation from the linear 
analysis (i Blaes fc Socrate^ l2001h . The Maxwell stress 
was found to be very similar to the value from sim- 
ulations with radiation pressure replaced by the same 
amount of gas pressure. Here we attempt to reproduce 
these results without invoking the FLD approximation. 
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as well as extend the study to higher ratios of radia- 
tion to gas pressure. The effects of the box size on the 
MRI turbulence will also be addressed fe.g.. iBodo et alj 
I2008D . In the case of gas pressure only, the properties 
of turbulence generated by the MRI have been studied 
extensively. For example, the ratio between Maxwell 
and Reynolds stress is always found to be ^ 4—5 
(jHawlev et al.lll995l:[Blacknian et al.ll2008D . and Maxwell 
stress is also fo und to scale with the magnetic pressure 
very well ( e.g., 'Blackman et al."2008; 'G uan et al.l 120091 : 
iHawlev et al. 2011; Sorathia ct al. 2012). Our goal is to 
investigate how these properties of MRI turbulence are 
affected by radiation, especially in a strongly radiation 
dominated flow. 

This paper is organized as follows. In Sj2l we describe 
the equations we solve and the shearing periodic bound- 
ary condition we use, especially for the radiation quan- 
tities. Then in ji3.1[ we repeat and ex tend the fiducial 
simulations done bv [Turner et al.l()2003f ). The saturation 
states of MRI for different total pressure are studied in 
Discussions and conclusions are given in Sj4l 

2. EQUATIONS 

We adopt the local shearing box approximation 
with shearing pe riodic boundary conditions (e.g., 
IHawlev et al.l[l995l ) for this work. In this approximation, 
we solve the equations of motion in a frame rotating with 
orbital frequency 51 at a fiducial radius from the central 
black hole. Curvature is neglected and the equations are 
written in local Cartesian coordinates (x, y, z) with unit 
vectors (i, j, k), which represent the radial, azimuthal 
and vertical direction respectively. Tidal and Coriolis 
forces are included in the local frame. 

The implementation of the local shearing box approxi- 
mation in Athena without radiation are described in de- 
tail by iStone fc Gardineii ()2010D . For radiation MHD, 
we need to add radiation source terms as described by 
JSD12. The radiation MHD equations in the mixed 
frame under the local shearing box approximation are 



dp 

dt 



V-ipv) = 0, 



djpv) 
dt 



+ V ■ipvv + P*) = -S,{P) 



-2pV,^qxi — 2Qk X pv, 



dE 
'dt 



V • [(S + P*)v - B{B ■ v)] = -cSr{E) 



fl^pv ■ {2qxi), 



dB 

'dt 



V X (v X B) = 0, 



dEr 
1 dFr 



V • Fr=cSr{E), 
-V-P, = Sr(P). 



C2 dt 

In the above equations, the shear parameter q = 
—dlnfl/dhir and g = 3/2 for Keplerian rotation while p 
is density, P* = {P + B'^/2)\ (with I the unit tensor), and 
the magnetic permeability p, = 1. The total gas energy 
density is 



(1) 



E^Eg + ^pv' + —, 



(2) 



where Eg is the internal gas energy density. We adopt 
the equation of state for an ideal gas with adiabatic index 
7 = 5/3. Then Eg = P/(7-l) and T = P/i?idcai/0, where 
Pidcai is the ideal gas constant. Radiation pressure, en- 
ergy density and fiux are P^, Er and Fr respectively 
while c is the speed of light. The radiation momentum 
and energy source terms are S r(P) and Sr{E). 

Following iJiang et al.l (|2012[ ) , we use a dimensionless 
set of equations and variables in the remainder of this 
work. We convert the above set of equations to the di- 
mensionless form by choosing fiducial units for velocity, 
density, temperature and pressure as ao, po, To and Pg^ 
respectively. A dimensionless parameter P can be de- 
fined as P = arT^/Pgfl, where is the radiation con- 
stant. Then units for radiation energy density E^ and 
fiux Fr are a^TQ and carT^. In other words, = 1 in 
our units. The dimensionless speed of light is C = c/oq. 
The original dimensional equations can then be written 
to the following dimensionless form 



d{pv) 

~dr 



| + V.(p.)=0, 
V • {pvv BB + P*) = -PSr(P) 



-2pil qxi — 2rtk x pv, 



dE 
'dt 



V • [(P + P*)v - B{B ■ v)] = -¥<CSr{E) 



Vi^pv ■ {2qxi), 



dB 

'dt 



V X {v X B)^0, 



dEr 

~W ' 

dFr 
"df 



:v-p,.=cs,(p), 
:v • p^==cSr(P), 



(3) 



where the source terms are. 



Sr(P) 



vEr 



+ ^(o'apP'* - CTasEr), 

Sr{E) = {(TapT^ - (JaEEr) 

, .V f vEr+V-Pr 
+ (CaF - CTsf) — • ll'r ^ 



(4) 



Frequency mean absorption and scattering opacities (at- 
tenuation coefficients) are daF, CTsf respectively while 
(T(jp and aaE are Planck mean and energy mean absorp- 
tion opacity (attenuation coefficients) . To change the di- 
mensionless equations to the dimensional form, we only 
need to replace C with c, set P = 1 and replace Fr with 
Fr/c. Note that we solve the time-dependent radiation 
momentum equations, and therefore include radiation in- 
ertia effect. 

The radiation pressure P^ is related to the radiation 
energy density Er with the variable Eddington tensor f : 



fP, 



(5) 



In principle, f should be calculated with our short char- 
acteristic module as described bv iDavis et al.l (pbl2i) . We 
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have tried this approach and find that the difference be- 
tween the Variable Eddington tensor (VET) and 1/31 is 
smaller than 10"'*. Thus for unstratified disk simulations 
with uniform density in the optically thick regime, the 
Eddington approximation is adequate, and we will adopt 
it in this paper. 

We use the recently developed Godunov radiation 
MHD method (JSD12) to solve those equations. This 
algorithm has bee n implemented and tested in Athena 
([Stone et al.l 1200 8^ as described in JSD12 with several 
improvements as explained in Appendix [X] The momen- 
tum and energy source terms resulting from the shear- 
ing box approximation a re added in the same way as in 
iStone fc GardineJ (|201Cll ). which are separated from the 
radiation source terms. 

2.1. Boundary Condition 

We adopt the usual shearing periodic boundary con- 
dition in X, and periodi c boundary conditio n in both 
y and z directions (e.g.. ^Haw lev et aD 11995^. Bound- 
ary conditions for the gas quantities are the same as in 
non-r adiative MHD simulations (e.g.. IStone fc Gardineii 
120101) . The radiation quantities Er and Fr are new vari- 
ables, and boundary conditions for them deserve detailed 
discussion. 

If the radial size of the simulation box is Lx, then the 
boundary condition for is 



Er{x,y,z) ^ Er{x ± L^, y T qVlL^t, z) 



(6) 



This boundary condition is the same as the boundary 
condition for density. Fr is the Eulerian radiation flux, 
which includes the co-moving and advected radiation 
flux. We first remap the co-moving radiation fiux and 
then add the advection flux at the new location. Thus 
the boundary condition for Fr is 

Fr^x i-> Fr^x{x ±Lx,yT q^Lxt, z), 
Fr,y i-> Fr^y{x ± Lx,yT q^Lxt, z) 
qVtLx 



T (1 + fvv) 



C 



Er{x±Lx,yT q^Lxt, z) 



Fr,z Fr^z{x ±Lx,yT q^Lxt, z), 



(7) 



where fyy is the yy component of the Eddington tensor 
at the new location. 

2.2. Energy Conservation 

Unlike the ZEUS code, Athena conserves energy to 
roundoff error when energies exchange between kinetic, 
internal and magnetic forms. However, as explained in 
JSD12, when photons and gas exchange energy, the ra- 
diation module in Athena does not conserve total energy 
to round off error because of our splitting scheme and 
the implicit matrix solver. With the special treatment of 
the energy source term described in JSD12, the energy 
error is kept small. This is particularly important for 
MRI simulations because the gas is continuously heated 
up during the simulations. If the energy error is not 
treated carefully, it can accumulate and eventually affect 
the dynamics. 

For unstratified simulations with shearing periodic 
boundary conditions, the change of the total energy in- 
side the simulatio n box is en tirely due to the work done 
on the walls (e.g., iHawlev et al.„1995i : .Gardiner fc Stonel 



l2005f ). With radiation, the volume averaged total energy 
is Et = {E + + PEr), where $ is the tidal potential 
$ = —qft^x^. The change of total energy is determined 

by 



LyL 



y^z J X 



{pVxSvy - BxBy) dydz = W{t), (8) 



dEt 
dt 

where 5vy = Vy + qflx, and the integral only needs to be 
done on one side of the radial boundary. The total work 
done on the simulation box within a certain time to is 
then Wt{to) = /o W(t)dt. To check the energy error, we 
calculate the energy change at to as 6E = Et{to) — Et{0), 
so that the energy error is SE — Wt{to). We calculate this 
energy error for all simulations and find that with respect 
to the total work Wt{to) within the first 100 orbits, the 
energy errors are always smaller than 1%. 

3. RESULTS 

In thi s section, we first re peat the fiducial simulations 
done by [Turner et al.l (j2003l ) and then extend them to a 
box with a larger radial size. In addition, we will study 
the saturation state of MRI turbulence with different ra- 
diation pressure. 

3.1. Simulations with Net Vertical Flux 



[Turner et all (|2003f) used the FLD module in ZEUS 
(|Turner fc Ston(:il2001f) to study MRI turbulence in the 
optically thick regime of a unstratified disk model. The 
radial size of the simulation box was fixed to be one scale 
height. Here we repeat their simulations with net vertical 
flux to compare the results from the two different codes. 
Because the whole simulation box is very optically thick, 
we expect to find similar results as those computed with 
FLD. We also extend the radial size of the simulation box 
to four scale height to study possible new effects with 
radia tion in this larger radial domain (e.g., IBodo et all 
[2003) . 

3.1.1. Initial Condition 

We take th e parameters for location I in table 1 of 
[Turner et al.l ([2003). This corresponds to a radius 67.8rG 
from a IO^Mq supermassive black hole, where the grav- 
itational radius rc = GM/c^ and M is the black hole 
mass. The whole simulation box has uniform density 
po = 2.89 X 10"^ g cm-3 and temperature Tq = 2.71 x 10^ 
K. The density and temperature correspond to the mid- 
plane in an a disk model with a = 0.01 and accretion rate 
10% of the Eddington rate. The mean molecular weight 
is assumed to be 0.6. Gas and radiation are in thermal 
equilibrium initially. We include both electron scatter- 
ing opacity kt = 0.4 cm^ g~^ and free-free absorption 

opacity = lO^V^M-Pg/lT " 1)]"^^^ cm^ g-i. Tak- 
ing the unit of density to be poi temperature to be Tq, 
pressure to be Pg,Q = po^ideaiTo (with i?idoai to be the 
ideal gas constant), time to be 1/51, and length to be 
Lq — {Pgfil pqY''^ /Vl, the dimensionless parameters in 
our equations are P = 376.3, C = 4895.6, (j^f = 1953.9, 
and aaF = <^aP = <^aE = 0.105. Thus the ratio between 
the background radiation and gas pressure is ~ 125. The 
ratio between radiation energy density and rest mass en- 
ergy of the fluid is arT^/ (poc^) = 3.14 x 10"^ Table 
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[T] lists other parameters for the five simulations in this 
section, with different box sizes and resolutions. 

The initial magnetic field is = By = 0, Bz = 0.0353. 
The most unstable MRI wavelength is = 2ttva,z/^ — 
2-K^fWJ~pl?L = 0.22Lo. Notice that with the unit of 
length we choose, the size of the simula tion box are 
smaller than the fiducial runs done by [Turner et al.l 
(p003i) . The net magnetic vertical flux we use is also 
smaller so that the ratio between the most unstable MRI 
wavelength and t he vertical size of t he simulation box 
is the same as in [Turner et all (pOOl . We have found 
that choosing a smaller field strength produces less vig- 
orous channel solutions (e.g.. iGoodman &: Xulll994 ) and 
leads to a more rapid transition to turbulence. Initial 
random perturbations are applied to both gas pressure 
and velocity field 



Pn 



S3 - qflx,Vz = S4, (9) 



where is a random number uniformly distributed be- 
tween —0.1 and 0.1 while 62,63,64^ are random num- 
bers uniformly distributed between —0.02 and 0.02. 
Athena's orbital advection scheme, as described in 
iStone fc GardineJ (|2010( ). is used for all the simulations 
to speed up the calculation. For the radiation field, the 
advective radiation fiux due to the mean orbital motion 
is separated from other terms and added as described in 
Appendix lA. 21 

3.1.2. Time Evolution 

The fiducial simulation in Tabic [U ZFRS, has the stan- 
dard geometry = Lo,Ly = 4Lo,Lz — Lq. To see the 
effect of the size of the simulation domain, we increase 
the radial size such that Lx = 4io a-nd label this simula- 
tion ZFRL. For comparison, the simulations ZFRNS and 
ZFNRL have the same parameters as ZFRS and ZFRL 
respectively, except that radiation pressure is replaced 
with the same amount of gas pressure. The differences 
between the two sets of simulations will be due to the 
radiation field. Simulation ZFRS64 has double the reso- 
lution of ZFRS. 

In Figure [Jl we show the time evolution of radiation 
and gas pressure, three components of magnetic and ki- 
netic energy densities as well as the total work done on 
the wall of the simulation box per unit time for the fidu- 
cial run ZFRS. Because photons cannot leave the sys- 
tem in the unstratified shearing box models, and more- 
over the box is continuously heated up due to the work 
done on the walls, there is no equilibrium solution and 
radiation and gas pressure continue to increase. Since 
Er oc T'* while gas internal energy density Eg (x T, radi- 
ation energy density increases faster and most of the work 
is converted to the radiation energy density. Although 
total energy cannot reach an equilibrium state, the mag- 
netic and kinetic energy densities saturate quickly after 
the initial transient. They fiuctuate around some mean 
values, which do not show any systematic change with 
pressure within the 100 orbits. This is also true for the 
Maxwell and Reynolds stress, the histories of which for 
four simulations are shown in Figure [2] Time histories of 
these quantities are very similar for simulations with or 
without radiation and it is independent of the box size. 

3.1.3. Statistical Properties 
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Fig. 1. — Time history of the spatially averaged quantities from 
simulation ZFRS with net vertical flux. The radiation (the red 
line) and gas (the black line) pressure are shown in the top panel. 
Components of the magnetic and kinetic energy densities are shown 
in the middle two panels. The black, red and blue lines represent 
X, y and z components respectively. Background shearing is not 
included in the calculation of kinetic energy density. The bottom 
panel shows the total work done on the wall per unit time, calcu- 
lated according to equation [8] 

In Table [2l we calculate different spatial and temporal 
averaged quantities for the five simulations. The tem- 
poral averages are taken between 10 and 100 orbits to 
avoid the initial transient before the MRI saturates. In 
that table, astross , ctmag and a are defined as 



Q^stress - 



ii-BxBy)) _ {{-BxBy)) 

((P«.<S«,))'"'"^«- {{By 2)) 
JJ^-BxBy^-{_pVxSv^ 
Po 



(10) 



where Pq is the initial total pressure and ((• )) represents 
spatial and temporal averaging. 

As shown in Table [2l the averaged Maxwell stress, 
Reynolds stress and magnetic pressure are very similar 
for simulations ZFRS and ZFRL, which suggests that 
changing the box size along radial direction does not af- 
fect the statistical properties of the turbulence. However, 
as shown in Figure [H when is increased by a fac- 
tor of 4, there are fewer large spikes compared with the 
fiducial run. This is also true for the pure gas pressure 
runs (sim ulations ZFNRS a nd ZFNRL). This is consis- 
tent with iBodo et al.1 ()2008f ) for MRI without radiation. 
The spikes in the fiducial run are due to the recurrent 
channel solutions. When the radial size is increased, 
the channel solution can be mo re easily destroyed ei- 
ther by the para sitic modes fe.g.. IGoodman fc Xulll99j : 
IPessah fc Goodman 200S), or b y the interactions be- 
tween different MRI modes (e.g.. iLatter et al.ll2009l ). 

When radiation pressure is replaced with the same 
amount of gas pressure. Maxwell stress for simulation 
ZFNRS is about 50% larger than for simulation ZFRS. 
The difference is primarily due to the large spikes in the 
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TABLE 1 

Parameters for simulations with net vertical flux 



Label L:^/Lo Grids/Lp Ts r„ Pg Pr -B^,o/2 



ZFRS 


1 


32 


1954 


0.11 


1 


125 


6.25 X IQ-"* 


ZFRL 


4 


32 


1954 


0.11 


1 


125 


6.25 X 10-4 


ZFNRS 


1 


32 






126 




6.25 X 10-1 


ZFNRL 


4 


32 






126 




6.25 X lO""' 


ZFRS64 


1 


64 


1954 


0.11 


1 


125 


6.25 X 10""' 



1 Ts is the electron scattering optical depth per Lq. 

^ Ta is the Planck mean free- free absorption optical depth per Lq. 



TABLE 2 

Temporal and Spatial Averaged Properties of Simulations with Net Vertical flux 



Label 


{{-B.By)) 




((5^2 » 


m/2)) 


Ctstress 


ZFRS 


0.154 


0.0282 


0.320 


0.0398 


5.46 


ZFRL 


0.162 


0.0291 


0.338 


0.0530 


5.57 


ZFNRS 


0.235 


0.0347 


0.518 


0.0560 


6.77 


ZFNRL 


0.123 


0.0273 


0.240 


0.0390 


4.51 


ZFRS64 


0.321 


0.0486 


0.761 


0.137 


6.71 


Label 


({B-il2)) 




((T/To » 


dmag 


a 


ZFRS 


0.268 


0.0113 


1.04 


0.481 


0.00145 


ZFRL 


0.263 


0.0220 


1.05 


0.478 


0.00151 


ZFNRS 


0.448 


0.0150 


1.61 


0.453 


0.00214 


ZFNRL 


0.188 


0.0130 


1.41 


0.513 


0.00119 


ZFRS64 


0.556 


0.0683 


1.09 


0.428 


0.00297 



A 
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V 




0.004 
0.002 


0.0025 
0.002 

0.0015 
0.001 

0.0005 








1 1 1 1 1 1 1 1 1 1 
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Fig. 2. — History of the spatially averaged Maxwell stress (loft panels) and Renolds stress (right panels) for the four simulations with net 
vertical flux. Parameters of the four simulations are given in Table [l] The stress is normalized with respect to the initial total pressure. 
The number in each panel is the time averaged value between 10 and 100 orbits. 



non-radiative run, which are weaker by almost a factor 
of 5 for the simulation with radiation. If we exclude 
those spikes, ZFNRS has a very similar saturation level 
as ZFRS. The change in temperature is 50% larger dur- 
ing the first 100 orbits for ZFNRS than ZFRS, although 
the total work done on the fluid is very similar for the 
two simulations. This is because the heat capacity for the 
photons is much larger than an ideal gas. When Lx is in- 
creased by a factor of 4, ZFNRL and ZFRL have roughly 



the same Maxwell stress a nd Reynolds s t ress. This is 
consistent with Figure 2 of [Turner et"an (p003l ). which 
shows that the stress is very similar when the net vertical 
flux and total pressure are the same for their simulations 
NV4 and RV4. The stresses f r om o ur simulations are 
smaller than what lTurner et ahl (j2003i) found because we 
use a smaller net vertical flux. This is roughly consistent 
with the scaling between the saturated Maxwell stress 
and the initial net vertical magnetic flux as found by 
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IHawlev et all (ITOOl. Simulations RV41i, RV4 and RV41 
in [Turner et al.l ( 20031 ) also have very similar stress, al- 
though the three simulations have quite different opaci- 
ties and thus different coupling between the photons and 
gas. 

The ratio between Maxwell and Reynolds stress, 
Q^stress, varies from 4.5 to 6.8 for the first four simulations 
in Table[2j This is a lso consistent with p revious MRI sim- 
ulations wi th (e . g.. [Turner et al.ll2003[ ) or without (e.g., 
iSano et al.ll2004l : iDavis et al.ll201 0^"radiation. When res- 
olution is doubled for simulation ZFRS64, the Maxwell 
stress is doubled while Reynolds stress is increased by 
70% compared with simulation ZFRS. But the ratio be- 
tween Maxwell stress and magnetic pressure amag varies 
in the small range 0.43 ~ 0.51 for all the five simulations. 
This is consistent with previous work, which found am ay, 
or th e tilt ang le defined based on it (e.g.. iGuan et al.l 
120091 IHawlev et al. 20fi|; iSorathia et aITl2012| ). to be 
roughly constant for converged simulat ions with differ- 
ent box sizes and initial co nditions (e.g.. lBlackman et al.l 
120081: IHawlev et alll2011[ ). Most previous MRI simula- 
tions without radiation find amag ~ 0.5 (e.g.. lGuan et al.l 
2009') ■ althou g h a r ange of 0.3 — 0.4 is also reported by 



Hawley et al.l ()2011[ ) for the stratified disk simulations 



3.1.4. Density fluctuation 

When there is no radiation field, turbulence gener- 
ated by the MRI i s almost incompressible. However, 
I Turner et al.l (|2003[ ) finds that if the photon diffusion 
length per orbit is larger than the most unstable MRI 
wavelength, the turbulence generated by MRI with ra- 
diation field becomes very compressible. This phe- 
nomenon can be characterized by the standard devi- 
ation of density 6p = ^ {{p — PoY). For simulations 
ZFRS, ZFRL and ZFRS64, photon diffusion length per 

1/2 

orbit [27rC/ (3(craF + f si?)ri)] ' = 2.3 while the most 
unstable MRI wavelength is only 0.25. Thus we are in 
the regime in which large density fluctuations are ex- 
pected. Indeed, as shown in Figure |3l for simulations 
without radiation, 5p is only 0.16%po for both ZFNRS 
and ZFNRL. While for the simulations with radiation, 
5p reaches 15%po for ZFRS and 11%pq for ZFRL. 

In summary, we find that in the very optically thick 
regime, our results are similar to those computed us- 
ing ZEUS w i th FL D for the parameters explored by 
I Turner et al.l (j2003[ ). Changing the box size does not af- 
fect the statistical mean values of the turbulence, which is 
true for simulations with or without radiation. Although 
the total pressure can increase by a factor of 2 within 100 
orbits, saturation levels of Maxwell and Reynolds stress 
do not show a systematic trend of change with pressure 
in the unstratified disk simulations. 

3.2. Simulations with Net Toroidal Flux 

In this section, we study the possible effects of radi- 
ation on the saturation states of the MRI for the cases 
with net toroidal flux. This conflguration of magnetic 
field is commonl y used for most stra tified accretion disk 
simulations (e.g.. IHirose "eralll2009l) . 

3.2.1. Initial Condition 

Here we carry out a set of simulations with different 
temperature as listed in Table |3l For these simulations. 
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Fig. 3. — Standard deviation of the density for the four simula- 
tions with net vertical fiux listed in Table [l] The number in each 
panel is the time averaged value. The first and the third panels are 
for simulations without radiation for comparison, which have much 
smaller density fluctuation compared with the other two panels. 

we choose our fiducial parameters to be the mid-plane 
va lues of the radiation dominated accretion disk studied 
bv IHirose etHI (|2009[ ). The disk is located at 30rG from 
a 6.62M0 black hole. The fiducial density and temper- 
ature are po = 5.66 x 10^^ g cm^^ and Tq = 2.63 x 10'' 
K. The orbital frequency Vl — 190 s~^ and the time unit 
is Mean molecular weight is assumed to be 0.6 and 

1/2 

the length unit is Lq = {Pg^/po) The dimension- 

less parameters P — 17.6 and C = 49 7 . We adopt the 
same opacities as used in IHirose et al.l (|2009[ ). Electron 
scattering opacity is 0.33 cm^ g~^, Planck and frequency 
mean free-free opacity are 3.7 x lO^^p^^^ [Pglil ~ 



and lO^V^/^ [PgKl 



-7/2 



gl\i -l)Y ■' cm' g " respec- 
the dimensionless electron scat- 



cm 

tively. In our units, 
tering and Planck mean absorption opacities (attenua- 
tion coefficients) with density po a-nd temperature Tq are 
= 5.94 X 10^ and o„p — 31.05. The box size is fixed 



to be Lx — Lq, Ly = 4Lo and — 4Lo- 

All the simulations listed in Table [3] start with a uni- 
form density 0.5po but we change the initial temperature 
to setup different radiation pressures. Initial perturba- 
tions are applied according to equation |9l For the fidu- 
cial parameters, the ratio between radiation energy den- 
sity and rest mass energy of the fluid is a^Tg/ (poc^) — 
7.13 X 10~^. For the case when temperature is 3ro, the 
ratio is increased to 5.78 x 10"'^. 

Simulations YFRl, YFR2 and YFR3 start from ther- 
mal equilibrium states with temperature To, 2To and 
3ro- The ratio between radiation pressure and gas pres- 
sure is 11.7, 93.9 and 317.0 for the three simulations re- 
spectively. Two oppositely twisted magnetic flux tubes 
with the same azimuthal flux are initially placed at 
— 0.35Lo < z < 0.35Lo. B^. and are initialized with 
the following vector potential Ab- 



ABix,y,z) 



-sign(z)^ [1 





cos (vrr)] 



r > 1 
r < 1 



7 



where 



{x/L.,f + {z- 0.35io)V(0-35Lo)2 
{x/L^f + (z + 0.35Lo)V(0.35Lo)^ 



1/2 



1/2 



Then By is initiahzed as By 



z > 
(12) 

z < 

1/2 



{Bl/2-Bl-Bl, 
such that magnetic pressure is uniform in the simula- 
tion box. We adopt -Bo/2 = 0.04 for aU the simula- 
tions in this section so that the net toroidal magnetic 
flux is always the same when we change the pressure. 
Simulation YFNR has the same setup as YFRl except 
that radiation pressure is replaced with the same amount 
of gas pressure. Simulations YFRlRes, YFR2Res and 
YFRSRcs arc doubled resolution version of YFRl, YFR2 
and YFR3 respectively. 




40 60 
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100 



Fig. 4. — The same as Figurc[T]but for simulation YFRSRes with 
net toroidal flux. 



3.2.2. Time Evolution 

Temporal evolution of radiation and gas pressure, three 
components magnetic and kinetic energy densities as well 
as the total work done on the wall of the simulation box 
per unit time for simulation YFR3Res are shown in Fig- 
ure m Because the initial radiation energy density is so 
large, within 100 orbits, the total energy only increases 
16.7% due to the work done on the wall of the simula- 
tion box. Similar to case with net vertical flux, kinetic 
and magnetic energy densities saturate quickly after ini- 
tial transient and fluctuate, although the total energy 
density never reaches an equilibrium value. 

In Figure [5l we show snapshots of p, E^.^ —B^By and 
{Frfi ■ 6v)/\dv\ for simulation YFRlRes at time 25.1 or- 
bits. Here Sv is the perturbed velocity with the back- 
ground shear subtracted. The same quantities for simu- 
lation YFR3Res at time 72.0 orbits are shown in Figure 
[HI Similar to the simulations in Table [I] YFRlRes and 
YFR3Res are also very compressible. The ratio between 
the maximum and minimum densities is ~ 10. Most of 
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Fig. 5. — Snapshots of simulation YFRlRes at time 25.1 orbits. 
From left to right, top to bottom, are density p, radiation energy 
density Er, Maxwell stress —BxBy and co- moving flux projected 
along the direction of perturbed velocity {Fr,o ■ Sv)/\&v\. 




Fig. 6. — The same as Figure [5] but for simulation YFRSRes at 
time 72.0 orbits. 



the Maxwell stress is located in regions where the den- 
sity variations are large, with the largest values in the 
low density regions. The radiation energy density has a 
similar distribution as density because Er is also com- 
pressed when p is compressed. Therefore the co-moving 
flux Frfi is also large when the density gradient is large. 
Radiation energy density E^ in simulation YFR3Res is 
81 times the value in simulations YFRlRes while the co- 
moving flux is almost a factor of 8 larger in simulation 
YFR3Res than YFRlRes. This means that the radia- 
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tion force on the fluid in YFRSRes is also about 8 times 
larger than the radiation force in YFRlRes. The sign 
of the quantity (.Fr.o • represents the relative 

direction between the co-moving flux and the velocity 
field. When it is positive, photons accelerate the fluid 
during expansion while it is negative, photons deceler- 
ate the fluid during compression. The net effect is that 
part of the kinetic energy is converted to gas internal and 
radia tion energy, resulting in damping due to radiation 
drag (lAgol fc Krolik "1998l The fact that co-moving flux 
in simulation YFRSRes is about 8 times larger than the 
co-moving flux in simulation YFRlRes, implies that the 
radiation damping effect in YFRSRes will also be much 
larger than this effect in YFRlRes. The consequences of 
the radiation drag effect will be examined in more detail 
in the following sections. 



TABLE 3 

Parameters for simulations with 
net toroidal flux 



Label 


Ps 


Pr 


N/Lo 


YFNR 


6.37 




32 


YFRl 


0.50 


5.87 


32 


YFR2 


1.00 


93.92 


32 


YFR3 


1.50 


475.47 


32 


YFRlRes 


0.50 


5.87 


64 


YFR2Res 


1.00 


93.92 


64 


YFRSRes 


1.50 


475.47 


64 



S.2.S. Statistical Properties 

Averaged properties between 10 and 100 orbits for 
these simulations are shown in Table SI By compar- 
ing YFRl, YFR2, YFRS and YFRlRes, YFR2Res and 
YFRSRes, we see that when the initial temperature is Tq, 
the simulation with S2 cells per Lq is already converged, 
while for the simulations with initial temperatures 2To 
and STq, this is not the case. To check the convergence 
of the high temperature runs, we carry out another simu- 
lation with 128 grids per Lq for initial temperature STq. 
We only run this high resolution module for 40 orbits. 
Averaging over the last 20 orbits, {{—BxBy)) = 0.220, 
{{pVxSvy)) — 0.00780 and amag — 0.258. Comparing 
this run and simulations YFRS and YFRSRes, we find 
that the differences between 64 cells and 128 cells per Lq 
are much smaller than the differences between S2 cells 
and 64 cells per Lq. For example, the Maxwell and 
Reynolds stress only change within the statistical fluc- 
tuations when we go from 64 to 128 cells per Lq. Thus, 
we infer that the high temperature runs with 64 cells 
per Lq are approaching convergence. The analysis be- 
low will focus on the four simulations YFNR, YFRlRes, 
YFR2Res and YFRSRes. 

S.2.4. Change of Maxwell stress with total pressure 

All the simulations done in this paper have fixed ver- 
tical box size and thus we cannot study the increase of 
Maxwell stress due to the increase of disk scale height, 
which we expect to be the primary source of connec- 
tion between stress and pressure. This will be stud- 
ied in future work. However, as Figure [7] shows, even 
with fixed box size, the Maxwell stress is increased with 
increasing total pressure. Compared with simulation 



YFRlRes with temperature T = Tq, initial total pres- 
sure is increased by a factor of 15 when temperature 
T = 2To (run YFR2) and a factor of 75 when tempera- 
ture T = STq (run YFRS). However, total stress is only 
increased by a factor of 2.67 from simulation YFRlRes to 
YFR2Res, and a factor of 5.18 from simulation YFRlRes 
to YFRSRes. 

To identify the cause of increase in Maxwell stress as 
total pressure increases, we have studied the magnitudes 
of the radiation source terms in detail. In equation 
the terms that are responsible for the radiation drag 
effect are the radiation momentum source term Sr{P) 
and the velocity dependent radiation work terms in the 
radiation energy source term Sr{E). The amplitude of 
these source terms for simulations YFRlRes, YFR2Res 
and YFRSRes can be measured directly from the sim- 
ulation data. We first calculate the spatial and tem- 
poral averaged radiation momentum source term Sp = 
{{F{aaF + crsF)\Frft\ ))• The temporal average is taken 
over 10 and 100 orbits to avoid the effect of initial tran- 
sients. For YFRlRes, YFR2Res and YFRSRes, in units 
of poaoil, Sp are 1.01, S.90 and 6.S8 respectively. As 
expected, the radiation momentum source term is in- 
creased when temperature and thus radiation pressure 
is increased. Second, we calculate tr = {{Ek ))/ {{Wr)), 
the averaged ratio between the kinetic energy Ek and 
the radiation work term Wr = —P{aaF — o'sf)'^ ■ Fr.o- 
This is the time scale to change the kinetic energy due 
to the work done by the radiation force. Here Ek does 
not include the kinetic energy due to background Kep- 
lerian motion. For simulations YFRlRes, YFR2Res and 
YFRSRes, in our time unit tr is —5.70, — S.SS and 
— S.20 respectively. Here U is negative, which means 
that the co-moving flux prefers to point to the opposite 
direction of fluid motion and decelerates it on average. 
The time scale is shorter for simulations YFR2Res and 
YFRSRes, only about half an orbital period. To make 
sure that the momentum source term Sr{P) is the key to 
the increase the Maxwell stress, we have also carried out 
another test in which we remove the momentum source 
term Sr(P) and the radiation work terms in the radi- 
ation energy source terms. We use the same setup as 
simulation YFRSRes. We find almost the same Maxwell 
stress and Reynolds stress as simulation YFRlRes and 
YFNR. Therefore, if there is no damping of fluid motion 
by the photons, we find the same saturation level as in 
the pure gas pressure case. This is strong evidence that it 
is radiation damping that is responsible for the increase 
of the Maxwell stress. 

The importance of the radiation drag effect c an also 
be qu antified by the photon bulk viscosity (e.g.. iCastoii 
I2004D . In o ur units, the r adiation bulk viscosity is (equa- 
tion 6.80 of lCastoil[200l 

^r = 4P£;,/ (9atC) , (IS) 

where at = (JaF + CsF • A Reynolds number can be de- 
fined based on the radiation bulk viscosity, and the typ- 
ical velocity, length scale and density (isothermal sound 
speed Cs, scale height Cs/0 and fiducial density po re- 
spectively), that is 

i?e = PQCI/ (n^r) . (14) 

For simulations YFRlRes, YFR2Res and YFRSRes, the 
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TABLE 4 

Temporal and Spatial Averaged Properties of Simulations with Net Azimuthal flux 



Label 


{(-B,:By)) 




((BV2» 


{{B'i/2)) 


*^strcss 


YFNR 


0.0330 


0.00730 


0.123 


0.0122 


4.52 


YFRl 


0.0469 


0.00580 


0.236 


0.0169 


8.09 


YFR2 


0.0737 


0.00290 


0.322 


0.0237 


25.4 


YFR3 


0.0547 


5.27 X 10-4 


0.261 


0.0129 


104 


YFRlRes 


0.0425 


0.00670 


0.199 


0.0175 


6.34 


YFR2Res 


0.125 


0.00660 


0.666 


0.0462 


18.9 


YFR3R0S 


0.249 


0.00570 


1.01 


0.0913 


43.7 


Label 






{{T/To » 


CKmag 


a 


YFNR 


0.107 


0.00410 


3.03 


0.267 


0.00632 


YFRl 


0.210 


0.00920 


1.21 


0.199 


0.00826 


YFR2 


0.281 


0.0175 


2.01 


0.160 


0.000806 


YFR3 


0.236 


0.0117 


3.02 


0.210 


0.000116 


YFRlRes 


0.174 


0.00790 


1.20 


0.214 


0.00772 


YFR2Res 


0.592 


0.0279 


2.09 


0.188 


0.00140 


YFRSRes 


0.861 


0.0561 


3.05 


0.247 


0.000534 



Reynolds numbers are 1.91 x 10^ 1.18 x lO'^ and 2.33 x 
10^ respectively. For simulation ZFRS, this number is 
5.71 X 10^. 

Our simulations also include grid scale dissipation set 
by numerical viscosity and resistivity. In order for the 
photons to control the magnetic dissipation rate, the 
Prandtl number defined by the ratio between radiation 
bulk viscosity and the numerical resistivity needs to 
be larger than 1 . Numerical resisti vity is likely prob- 
lem dependent. iSimon et al.l (j2009r ) quantifies the ef- 
fective magnetic Reynolds number based on numerical 
resistivity in Athena t o be around 5 x 10^ — 1 x 10''. 
ISimon fc HawlevI ()2009f ) also find that when the Reynolds 
number for the explicit (shear) viscosity is smaller than 
6400, the Maxwell stress begins to show a significant in- 
crease. Therefore we suggest that radiation pressure in 
simulations YFR2Res and YFRSRes is large enough to 
change the Maxwell stress significantly while in simula- 
tions YFRlRes and ZFRS, it is not. 

A similar result for the unstratified disk model was also 
report ed for the pure gas pressure case by Sa no et aD 
(j2004D . They studied the saturation level of MRI for 
a wide range of gas pressure and found that Maxwell 
stress has a weak dependence on the gas pressure. 
For the zero net flux case, they found {{—BxBy)) oc 
P^/'^ while for the net vertical flux case, they found 
{{-B^By)) cx Pi/s. There is no definite explanation 
for the weak dependence of Maxwell stress on the gas 
pressure. They argue that higher gas pressure can re- 
duce the magnetic reconnection rate and thus reduce 
the dissipation rate of magnetic energy, resulting in a 
higher Maxwell stress in the saturation state. More re- 
cently, studies with explici t viscosity and resistivity (e.g. , 
Lesur fc Longaretti 20 071: iFromang fc Papaloizoul 120071 : 
Simon fc HawleVii2009i : ISimon et al.ll2009( ) foundlhat at 
moderate Reynolds number Re, the Maxwell stress from 
the saturation states of MRI turbulence increases with 
the magnetic Prandtl number Pm = v/r], which is the ra- 
tio between microscop ic visc osity and resistivit y. Follow- 
ing Eilbus^H^Ie^ (Hill) , ISimon fc HawlevI (pOOl ) ar- 
gues that viscosity can prevent motion that would bring 
field together on small scales and therefore reduce the 
magnetic reconnection rate. The increase of Maxwell 
stress with radiation pressure found here, may be due to 



a similar mechanism, except that the change of effective 
Prandtl number is caused by the damping of radiation 
field on the compressible motions in the fluid, rather than 
explicit viscosity. 

3.2.5. The ratio between Maxwell stress and Reynolds 

stress 

The change in the effective Prandtl number due to ra- 
diation drag also seems to affect the ratio of the Maxwell 
to Reynolds stress. In Figure III when the Maxwell 
stress is increased with increasing radiation pressure, the 
Reynolds stress stays at almost the same level. There- 
fore, the ratio between Maxwell stress and Reynolds 
stress is increased when radiation pressure is increased. 
As calculated in Table |3l this ratio is 6.34 for YFRlRes 
and it increases to 18.9 and 43.7 for YFR2 and YFR3 
respectively. As a comparison, the ratio is only 4.52 
for the simulation YFNR without radiation. This ra- 
tio reported in most MRI simulati ons with or without 
radiation is also around 4 — 5 (e.g., 'Hawlev et al."1995t 
iTurner et a l. 2003; San o et al.ll^ 04; Davis et al. 2010). 
However, when Si mon fc Hawlevr(|2009 ') includes explicit 
viscosity in their MRI simulations, they finds that the ra- 
tio between Maxwell stress and Reynolds stress increases 
when the stress and viscosity is increased, although the 
change of this ratio is not as large as what we find here. 

In the test run in which we remove the radiation mo- 
mentum source term described in the last section, we 
find the usual ratio between Maxwell stress and Reynolds 
stress as in the pure gas pressure case. We also notice 
that for all the simulations listed in Table [U where we 
find the radiation damping is not important, the ratio 
between Maxwell stress and Reynolds stress is also nor- 
mal. This ratio is only significantly increased for those 
simulations listed in Table El where the Maxwell stress 
is much larger than what the values would be if the ra- 
diation pressure is replaced with the same amount of 
gas pressure. Those results suggest that the change of 
the ratio between Maxwell stress and Reynolds stress is 
due to the same physical mechanism responsible for the 
increase of Maxwell stress discussed in the last section, 
which is the damping of the compressive motions due to 
radiation drag. We have also checked that this result is 
independent of field geometry. When we use net vertical 
magnetic fiux, we get the same result. 
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We also notice that a,nag for those simulations listed 
in Table |3] is systematically smaller than amag from non- 
radiative MRI simulations. It is also smaller than a,nag 
for the simulations listed in Table [l] which have net ver- 
tical magnetic flux. For simulations YFRlRes, YFR2Res 
and YFRSRes, amag varies from 0.188 to 0.247 with the 
average value to be 0.21. For simulations ZFRS, ZFRL, 
ZFNRS, ZFNRL and ZFRS64, a^ag is always between 
0.4 and 0.5, which is also the usual valu e reported by 
most non-radiative MRI simulations fe.g.. ^Hawle v et all 
120 111 ). A smaller a,nag is also found from stratified sim- 
ulations of radiation dominated accretion disk done with 
ZEUS (Krolik et al., 2012, private communication). 



3.2.6. Radiative Viscosity 

In the radiation dominated flow around the black holes, 
it has been argued that radiative viscosity can be impor- 
tant to transfer the angular momentum, which can even 
be the domi nant mechan ism for spherical accretion flows 
(e.z.. iLoeb fc LaoH 119921 ). Radiative viscosity will only 
be important when the radiation energy density is sig- 
nificant. Here, we can compare the viscous stress due to 
radiation with the Maxwell and Reynolds stress for simu- 
lation YFRSRes, which has the largest radiative viscosity 
in our simulations. 

In principle, the off-diagonal component of the co- 
moving radiation pressure in the x — y plane, which is 
responsible for the transfer of angular momentum, comes 
from two parts in our formula. The first part is the off- 
diagonal component of the Eddington tensor in the Eu- 
lerian frame. We have calculated the VET for simulation 
YFRSRes with the short-characteristic module, and find 
that the off-diagonal components of the VET is smaller 
than 10~*. Thus it can be neglected. The second part 
is the radiative shear viscosity, the formula of which is 
(e.g.. [Castor. 2004) 



IScTfC 



dvx 
dy 



dvy 
dx 



(15) 



This is caused by the velocity dependent radiation mo- 
mentum source terms. For simulation YFRSRes, the 
spatially and timely averaged value is {{Pr,xy )) = 4.13 x 
10"** in our fiducial units. The radiative viscosity is dom- 
inated by the background shearing. Therefore, even for 
simulation YFRSRes, the off-diagonal component of the 
co-moving radiation pressure is only 1.66 x 10~^ of the 
Maxwell stress and 7.2% of the Reynolds stress. 

4. DISCUSSION AND CONCLUSION 

With our recently developed Godunov ra diation MHD 
code, we have confirmed the results of iTurner et al.l 
(|200S[ ) for unstratified optically thick mid-plane of a ra- 
diation dominated accretion disk. In particular, we find 
when the photon diffusion length per orbit is larger than 
the most unstable MRI wavelength, the MRI turbulence 
becomes very compres sible. Furt hermore, we extend the 
fiducial simulations of ITurner et a l. (2003) by increasing 
the radial size by a factor of 4, which allows more radial 
modes to exist in the simulation box. The saturation 
level of the MRI from the larger simulation domains is 
very similar to what we find in the fiducial runs. Recur- 
rent channel solutions observed in small box simulations 
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Fig. 7. — History of the spatially averaged Maxwell and Reynolds 
stress for simulations with net azimuthal flux. Maxwell stress is 
plotted as the upper thick lines while the Reynolds stress is plotted 
as the lower thin lines. Parameters of all the simulations shown in 
this figure are listed in TableE] YFRlRes, YFR2Res and YFRSRes 
differ in the initial ratio between radiation pressure and gas pres- 
sure, which increases by a factor of 27 from YFRlRes to YFR,3Res. 
YFNR is a comparison simulation without radiation field but the 
same total pressure as YFRlRes. 



do not exist when radial size i s incr eased. This is con- 
sistent with what iBodo et al.l (|2008l ) found without ra- 
diation. In this sense, radiation pressure plays the same 
role as gas pressure in the optically thick regime. 

When the radiation energy density is increased to a 
non-negligible fraction (> 1%) of fiuid rest mass energy, 
we find that Maxwell stress is increased with increasing 
radiation pressure. At the same time, Reynolds stress 
seems to be independent of radiation pressure. There- 
fore, the ratio between Maxwell and Reynolds stress is 
increased with increasing radiation pressure. We suggest 
that this is due to the damping effect of radiation drag 
which introduces an effective bulk viscosity ^r- To quan- 
tify the importance of the radiation damping effect, we 
have calculated the Reynolds number based on (equa- 
t ion ITS]) and find values of ^ 2.33 x 10"^ for the simulation 
with the largest radiation pressure. This number de- 
pends on VEr/C in our units for a fixed opacity, which is 
actually the momentum of the photons. This can explain 
why we do not see the change of Maxwell stress as well 
as the ratio between Maxwell and Reynolds stress com- 
pared to the pure gas pressure case for the simulations 
listed in Table [H Although the ratio between radiation 
pressure and gas pressure is 125 for the simulations listed 
in Table m ¥Er/C is only 7.7%. While for the simula- 
tions listed in Table H FEr/C is 0.035, 0.567 and 2.87 
for T = To,2To, STq respectively. 

Because there is no cooling in the unstratified simu- 
lations, the disk is continuously heated up due to the 
dissipation of MRI. However, as the histories of these 
simulations show (Figure [5] and Figure [7]) , the secular 
build-up of energy in these simulations does not affect 
the diagnostics we are interested, such as Maxwell stress 
and Reynolds stress. This is not surprise, because the 
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effects discussed this in paper requires T to increase by 
a factor of 2 and 3, while during the 100 orbits of YFRl, 
YFR2Res and YFRSRes, T is only increased by 40%, 9% 
and 4% respectively. 

Even when radiation energy density reaches 1% of the 
fluid rest mass energy density as in simulation YFRSRes, 
we find that the off-diagonal component of the co-moving 
radiation pressure is much smaller than the Maxwell and 
Reynolds stress. This confirms that transfer of angular 
momentum by radiative viscosity is not the dominant 
mechanism when the turbulence due to MRI exists. 

We perform two sets of simulation with both high 
[N/Lq = 64) and low {N/Lq = 32) resolution. We find 
consistent (statistically equivalent) estimates for turbu- 
lent properties, such as the stress, which suggests these 
simulations are converged. Since we do not include ex- 
plicit resistivity, our results may have been sensitive to 
the effective amount of "numerical r esistivity" , which i s 
strongly resolution dependent (e.g.. iSimon et al.ll2009[ ). 
Our result implies that effective resistivity is sufficiently 
small that turbulent properties are set primarily by dy- 
namics at larger scales and not strongly impacted by the 
precise properties of grid scale dissipation. 

In all the simulations, we have fixed the the Eddington 
tensor to be l/3Sij in the Eulerian frame. Since our sim- 
ulations are approximating the mid-plane of an optically 
thick accretion disk, the assumption of isotropy is well 
justified, but it is generally expected that the radiation 
field will be nearly isotropic in the co-moving frame. As 
the difference between co-moving radiation pressure and 
Eulerian radiation p ressure in our units is v-Fr^/C (e.g., 
iMihalas fc MihalasI IT983) . the difference between Eule- 
rian Eddington tensor and co-moving Eddington tensor 
will be V ■ Fr.o / (CEr), which we confirm is less than 10^^ 
for all the simulations presented here. Therefore, this as- 
sumption should not appreciably impact our results. A 
test with Eddington tensor to be l/3Sij in the co-moving 
frame does confirm this conclusion. 

We have also calculated VET with the short charac- 
teristic module described in iDavis et al.l (|2012[ ) for the 
simulations shown in Table [T] For the unstratified shear- 
ing box disk models, there is no preferred direction for 
the specific intensity because periodic or shearing peri- 
odic boundary conditions are applied for all three di- 
rections. Furthermore, photon mean free paths for all 
the simulations are smaller than the cell size. Therefore, 
VET returned by the short characteristic module only 
differs from 1/31 by 10"'', when 24 angles are used. This 
confirms that Eddington approximation is valid for the 
unstratified disk models, which are designed to study the 
optically thick mid-plane of accretion disks. 

The increase of Maxwell stress with increasing radia- 
tion pressure and the change of astrcss have important 
implications. This study confirms that the damping ef- 
fect of large radiation pressure will not decrease the to- 
tal stress. Instead, the total stress is actually increased 
slowly with total pressure when radiation energy density 
is significant. 

These results also confirm the important role radiation 
viscosity plays as a source of microscopic dissipation in 
radiation dominated regions of accretion fiows. Although 
radiation viscosity always plays a subdominant role to 
Maxwell and Reynolds stresses in angular momentum 
transport, it will generally exceed the standard plasma 



viscosity. In cgs units, Eciuation (jl3p corresponds to a 
kinematic viscosity with 

An T** 
^i^= 3.4 X 10-25 i^cm^s-i. (16) 
9ktp c 

This can be comp ared with Equation (2) of 
iBalbus fc Henril (|2008[ ). which gives an estimate of 
the plasma viscosity Vp oc T^l"^ j p. In the radiation dom- 
inated regimes of accretion flows, > Vp when T < 10® 
K and it is often the case that Vr ^ Vp. It is still true 
for the hottest radiation dominated flows (T > 10® K) 

when p < 4.74 x 10'' (T/IO® K)^^^ g cm-^. So radiation 
is the dominant viscosity associated with microscopic 
dissipation. This makes radiation dominated accretion 
disk one of the few regimes where the physically relevant 
value of the viscosity can be numerically resolved with 
current computing resources. The ratio Vr/vp is larger 
for the accretion disks around supermassive black holes 
than solar mass black holes, which may lead to different 
properties of the accretion fiows for systems in the 
two different scales if the radiation viscosity becomes 
important. 

Despite the dominance of radiative viscosity, it is gen- 
erally true that plasma resistivity dominates over radia- 
tive resistivity, even in radiation dominated accretion 
flows IA koI & Krolik 1998). Therefore, Equation (1) of 
IBalbus fc Henril (|2008f l still gives the relevant magnetic 
resistivity rj. Since this quantity is generally very small 
compared with i^r, the corresponding magnetic Prandtl 
number Pm = fr/?? 3> 1 in radiation dominated flows. 
This conclusion is interesting in light of speculation that 
Pm may be the primary parameter that controls angular 
momentum transport by the MRI (Fromang et al. 2007). 
A large Pm would place radiation dominated accretion 
flows flrmly in the efficient angular momentum regime if 
this supposition is correct. 

In all the unstratified disk simulations, the total energy 
can never reach an equilibrium value because there is no 
cooling in these simulations. However, the Maxwell and 
Reynolds stress saturate quickly after an initial transient. 
Even when total pressure is increased by a factor of 2, 
the stress does not show any systematic trend of change 
within the first 100 orbits of the simulations. This is be- 
cause the saturation level in the unstratified disk models 
is lirnited by the box size. As found by Hawley et al] 
((T995I) . Maxwell stress in the saturated state increases 
linearly with box size in their unstratified disk models. 
Therefore, the unstratified disk models cannot be used to 
check the assumed relation between stress and pressure 
in the a disk model, which requires the change of disk 
scale height to connect them. The fact that the effec- 
tive a value decreases when total pressure increases for 
these unstratified simulations is an artifact of the closed- 
box assumption, which should not be interpreted liter- 
ally. When the total pressure increases, the disk scale 
height will also be increased, which gives more space 
for the turbulence to develop. The weak dependence 
of Maxwell stress on the radiation pressure due to the 
mechanism suggested here is probably not the dominant 
process in real accretion disks. More realistic stratified 
disk models will be studied in the future. Nevertheless, 
the phenomenon we find here is also likely to operate in 
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the mid-plane of the accretion disk near the Eddington 
limit. In the future, we will perform more realistic strat- 
ified disk models that will be better suited for addressing 
these questions. 
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APPENDIX 

IMPROVEMENTS TO THE NUMERICAL ALGORITHM 

The code used in this paper is based on the algorithm described in JSD12, with several important improvements. 
Based on our tests, those improvements are necessary to reduce numerical diffusion in the optically thick regime while 
leave the optically thin regime unaffected. 



Effective photon propagation speed 

The first change is the HLLE Riemann solver for the radiation subsystem (the last two lines of equation . As 
described in §3.3.2 of JSD12, the HLLE Riemann solver is used to calculate the flux for the radiation subsystem. In 
JSD12, the characteristic speed along each dire ction is always choseii to be \/JC, where / is the component of the 
E ddington tensor alon g that direction, following ISekora fc Stond (j201C0 (see their equation 46). However, as discussed 
by I Audit et al.l (j2002f) . when the optical depth per cell is larger than one, the numerical diffusive flux due to the HLLE 
Riemann solver can be much larger than the phy sical radiative flux, which can cause unphysical solutions. This issue 
does not appear in the tests done by iSekora fc Stona ()201Ql because they only consider pure absorption opacity. In 
the optically thick regime, the evolution of radiation energy density is dominated by the source term Caa{T'^ — Er) 
because the thermalization time is usually very short. However, for scattering opacity dominated cases, the energy 
source term is very small and the radiation flux term V • controls the change of Er . If the numerical diffusive flux 
is much larger than the physical radiation flux, the photon diffusion time will be signiflcantly underestimated. 

lAudit et all (^2002) fixes this issue by changing the term V • Pr in their dimensionless units (see their equation 46). 
This is equivalent to reducing the characteristic speed of the radiation modes in the optical thick regime. To find the 
appropriate effective speed for different frequency mean opacity at = <JaF + (^sF , we carry out a linear analysis for the 
following ID equations 

^ + /C^ = -Ca,F.. (Al) 
ot ox 

The background state is a static medium with a uniform radiation field. For mode with wavenumber fc, the propagation 
speed is 



r2 



^e// = V7cyi-^. (A2) 

In the optically thin limit when at 0, the effective speed of light Cg// approaches ^/JC In the optically thick limit, 
Ceff is much smaller than v7C. Especially, for modes with at /{4,fk'^) > 1, the propagation speed of the photons is 
zero, which means photons only diffuse through the medium instead of free-streaming. Because we solve the radiation 
subsystem implicitly, the distance that photons can travel for a give time step dt is not limited by the size of each cell. 
The lower limit for the wave number k should be l/(AiC). However, this value is usually too restrictive and makes 
the resulting matrix from the implicit radiation subsystem very hard to invert. Instead, another characteristic length 
is the size of each cell. Based on our tests, we find that if we calculate the effective propagation speed with a wave 
length of ten cells, the numerical diffusion is small enough not to affect the solution. If the size of each cell is Al, to 
connect the optically thin and optically thick regimes smoothly, the effective propagation speed we use is 



C:// = C^^ [l_exp(-T)], 

T={10Alxat)y{2f). (A3) 



In the optically thin limit, the effective speed is reduced to equation I A2 1 with wavenumber k — 1/(10A/). In optically 
thick limit, the effective speed C^jf oc 1/(AZ x at). The purpose of using C*jj instead of C is to reduce the numerical 
diffusion so that we can capture the physical radiation flux, which is small in the diffusion regime. As long as physical 
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radiation flux is dominant over the numerical flux in all regimes, the solution is more accurate. Tests presented below 
do confirm this. 

Advection Flux 

In the optically thick regime, the diffusive flux is very small, and the change of radiation energy density can be 
dominated by the advection flux. We find that the first order treatment for the radiation subsystem described in 
JSD12 is too diffusive for the advection flux. In order to improve the accuracy, we separate the advection flux from 
the diffusive flux and treat the advection flux explicitly with second order accuracy. 

For a given velocity v, the change of Er due to the advection flux is 



ot 



(A4) 



Because the radiation subsystem is separated from the MHD part, the flow velocity v is taken to be constant during 
this step. We first calculate the flow velocity at the edge of each cell Vi+i/2 by taking the average of cell centered 
velocities in the two neighboring cells. We use the van Leer slope (equation 48 and 49 of iStone fc Normanl[T99l) 
to interpolate E,. from cell centers to cell edges £'r,i+i/2- Then the advection flux on each cell edge is just 

= ^+l/2-E'r,z+l/2- 



The New Implicit Backward- Euler Step 

In order to incorporate the above two improvements into our implicit Backward Euler step, we first rewrite the 
radiation subsystem into the following form 



dEr 



■CV 



vEr 



V ■ Pr 



c 



V ■iv-Pr) = -V ■ (vEr) + CSriE), 
-CV - P^=CSr(P). 



(A5) 



Notice that Frfl = Fr — {vEr + v ■ P^.) /C is actually the co-moving flux. In this way, the advection part V • 
{v ■ Pr + vEr) is separated from the fl ux in the co- moving frame V • i^r,o- The advection part V • (t> • P r + v Er) is 
the advective radiation enthalpy (e.g.. lCastoHl200^ . Only the V • [vEr] part is treated as described in ^A.2\ while 
the V • (w • Pr) part is solved implicitly with the source terms. 
First, with the velocity fleld v and E^ at time step n, we calculate the advection flux along each direction J-i+i/2,j.ki 

^i,j+i/2.ki ^i.j,k+i/2 described in §A.2l for each cell. Then we calculate the implicit HLLE fluxes F^'"'^^ according 
to equation (39) of iSekora &: Ston"^ ([2010 ) for the following equation 

dEr / „ vEr + v - Pr 

~dt 



CV 



= 0, 



dFr 



:v • P. = 0. 



(A6) 



The following changes are made during the calculation of the HLLE fluxes. The maximum and minimum wave speeds 
we actually use are C*j^ and —C*jj respectively, where C*j^ is calculated according to equation I A3I In order to make 
the implicit matrix easier to invert for the cases with sharp opacity jump, the opacity at we use to calculate C*^^ on 
the cell edges is the average opacity of two neighboring cells. The term V • (v • P^) is calculated implicitly as 

V-{v- Pr) (i, J, k)=[v P^+i(* + l)-v P^+i(z - 1)] /(2Ax) 
+ [v ■ P^+i(j + l)-v P^+i(j - 1)] /(2Ay) 

+ [v ■ P^+i(fc + l)-v- P^+i(fc - 1)] /(2Az). (A7) 
Then radiation energy density Er is updated as 



e: 



r,i.j,k 



=e: 



At 
Ax 

HLLE 



■ t,j,k+l/2 



r.ij^k 

At 

At 
A^ 

AtV ■ {V-Pr 



TT^HLLE 
-■^ i+l/2j",fe 



TT'HLLE 

^ 2-1/2, 



At_ r. 
Ay 



T7.HLLE 

ij + l/2,fe 



TT'HLLE 
-■^ jj"-l/2,fc 



^HLLE 
■ i,j,k~l/2 



At . 

A^ L-^+i/2...'c 



[•^ ij+i/2,k ^ J'i,j-l/2,k] 



Sr{E). 



At 
'a'z 



k+l/2 



jj,fe-l/2j 



(A8) 



The special treatment of the energy source term Sr{E) described in JSD12 is also used here in order to reduce the 
energy error. The implicit equation for radiation flux Fr is unchanged. 
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Numerical Tests of the Improved Algorithm 

We find that the most useful tests to demonstrate the necessity of those improvements is the dyna mic diffusion test 
with p ure scattering opacity. The improvements do not change the results of other tests done by iSekora &: Stonel 
(|2010D and JSD12, which will not be repeated here. With pure scattering opacity, the radiation energy source term is 
zero and the change of radiation energy density is totally controlled by the divergence of radiation flux term. Inaccurate 
treatment of this term will easily show up in this case. In contrast to the static diffusion case, the dynamic diffusion 
test also requires accurate treatment of the advection term. 




Fi g. 8. — Dynamic diffusion test for the radiation subsystem with our improved algorithm. From left to right, the parameters for equation 
iMlare: at = 40000,!^^ = 40; at = 40, i/^ = 0.002; at = = 0.002. The initial gaussian profiles are all centered at x = 0. The black 

dots not centered at x = 0, the red dots and the green dots are solutions from earlier time to later times with time intervals shown in each 
panel. The blue lines are analytical solutions (equation I A9| | at corresponding times, which agree with the analytical solution very well for 
the regions not close to the boundaries. 

In the optically thick regime, the radiation subsystem can be approximated by diffusion equations. In ID with a 
constant opacity at and the Eddington approximation, the diffus ion equations and the analytic solution with initial 
Gaussian pulse for a static medium are given by equations (101) of lSekora &: Stond (|2010l ). except that at only contains 
scattering opacity here. For the dynamic diffusion test, we let the background fluid move with a constant velocity v. 
In the co-moving frame, the diffusion solution is the same as in the static diffusion case. In the Eulerian frame, the 
solution is 

with the mitial condition 

i;,(x,0)=exp(-i.V), Fr{x,0) = ^Er{x,0) + —Er{x,0). (AlO) 

Here D = C/{3at) is the diffusion coefficient and v is a. free parameter to control the shape of initial Gaussian pulse. 
We can initialize the Gaussian pulse in the code and compare the analytical and numerical solutions at time t. We 
use periodic boundary conditions so that we do not need very large simulation box when the pulse is advected. The 
following parameters are chosen for all the tests: C = 514.4, v = 1. We test the new algorithm in three cases with 
different opacities at = 40000,40,1. For the first case, the pulse is located within < 0.5 initially while for the 
other two cases, the pulse is located within \x\ < 50 initially. Results are shown in Figure [SI which confirms that the 
improved algorithm can reproduce the diffusion behavior accurately. We have also tried that if we do not use those 
improvements, the diffusion time is too short compared with what it should be for the pure scattering case. 
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